Non Standard Finite Difference scheme for Mutualistic Interaction 

Description 

Gianluca Gabbriellini 
January 4, 2012 

Abstract 

One of the more interesting themes of the mathematical ecology is the description of the mutualistic 
interaction between two interacting species. Based on continuous-time model developed by Holland and 
DeAngelis 2009 for consumer-resource mutualism description, this work deals with the application of 
the Mickens Non Standard Finite Difference method to transform the continuous-time scheme into a 
discrete-time one. It has been proved that the Mickens scheme is dynamically consistent with the original 
one regardless of the step-sizes used in numerical simulations, in opposition of the forward Euler method 
that shows numerical instabilities when the step size overcomes a critical value. 

1 Introduction 

The mutualism is acknowledged as a fundamental process in ecological systems: it is an interaction between 
individuals of different species, from which the interacting populations have beneficial effects. Like predation 
and competition, is now recognized as a consumer-resource (C-R) interaction [1]. The C-R interaction relates 
the process of nutrient transfer between a consumer and a resource: is a mechanism that describes the way by 
which individuals interact with each other [::]. All the types of mutualism belong to one of the three general 
forms of C-R interaction: uni-directional, bi-directional and indirect C-R mutualisms [3] . 

From a mathematical point of view only few authors have been attempts to modelize the mutualism, 
probably because it is difficult to retrieve the experimental data. One of the first attempts is based on 
Lotka-Volterra competition equations [s] , but it conducts to the boundless growth of the interacting population. 
Subsequently, other models are proposed where a stable equilibrium for populations are reached [•")], [(>]. 
Holland and DeAngelis [2] developed a general continuous-time model to describe the population dynamic for 
two mutualistic species; whose model is based on a modified version of the Rosenzweig-MacArthur [ ] model 
of predation. 

Although a continuous-time models based on differential equations constitute the first step to modelize 
a population evolution problem, for the more complexes issues investigated in mathematical biology, the 
numerical method is the only practicable way. Also, discrete-time models become essential when one wants 
to describe experimental data that have been collected with a certain interval of time. One of the critical 
aspects of the discretization methods is just the dipendence from this interval, the time step. Since the 
time step should be selected only in relation to the characteristics of the problem under examination, it is 
necessary to choose a reliable discretization method that allows to make this transformation without any 
restriction on time step, as well as not to introduce artifacts, change the linear stability properties etc.. Since, 
in general, these requirements are not satisfied from all the discretization methods (see i.e. [15]), in this paper 
is discussed the application of the Non Standard Finite Difference scheme [12] to the system of differential 
equations discussed by Holland and DeAngelis in the articles [■">] and [1]. Of this work is only considered 
the equations describing the population dynamic for bi-directional C-R mutualism, since the uni-directional 
problem is easily solved once it has solved the bi-directional one. 



1 



2 Time continuos approach 




A mathematical model for population growth can be generally described by a system of first order of 
autonomous ordinary differential equations: 

(1) 

with / : M" — > M" is supposed to be sufficiently smoothed; ^ — (,{t) : [0, +oo) — > M" are the coordinates, 
with initial condition € I^"; finally k = (fci,fc2,...) represents the system parameters. If / is globally 
Lipschitz then there exists a unique solution ^ for all t > tg, hence the eq. (1) defines a dynamical system on 
M". For a system defined with eq. (f ) we define a steady state solution (or also fixed point of /) a point 
f £ M" that satisfies the relation: 

/(I) = (2) 

In the continuation of the text will be used the set Fc = {^|/(^) = 0, ^ G M"} to refer to the set of steady 
states solutions. 

If / € it is possible to calculate the n x n Jacobian matrix as here indicated with Jf . In term of the 
Jacobian matrix, as given (j{Jf) the set of its eigenvalues, a theorem ensures that a steady-state ^ € M" for 
which J has no eigenvalues with zero real parts (hyperbolic steady-state) is [ , ] : 

• asimptotically stable if and only if 3ffA < 0, for all A G a[Jf), 

• unstable if and only if 5RA > 0, for all A G a{Jf). 
The same is not true for a non-hyperbolic steady-state. 



3 Finite difference approaches 

To transform a continuos-time model into a discrete one, the continuos variable t G [0, oo) must be replaced by 
the discrete variable n G N and the variable ^ must take discrete values The result is a difference equation. 
Let / : M" — > M" consider a sequence {^nj^o- can be defined by a mapping A : M" x M" — > M" of 
the form i?(^„+i,^„). In some case, such as the one discussed on this paper, is possible that ^n+i is given 
explicitly in terms of 

C„+i = F{i^) (3) 

where F : C M" — ^ M". 

For a dynamical discrete system defined by eq. (3) a steady-state ^„ G M" respects the following 
conditions: 

Filn) = L (4) 

Likewise to the continuous case, also in discrete case it is useful to define with r^; = {^n\P{^n) 0, G 11^"} 
the set of steady-states. 

To check the stability of the system described by eq. (4) the procedure is similar to the continuous case. 
If G is possible to calculate the n x n Jacobian matrix, indicated with Jp- Given (j{Jf) the set of the 
Jacobian eigenvalues, a theorem (see i.e. [ ]) ensures that a steady-state ^„ G K" is: 

• locally asymptotically stable if and only if 5RA < f VA G (t{Jf)'- this point is an attractor. 

• unstable if and only if JRA > I VA G a{JF)'- this point is a repeller. 

• no conclusions on stability if 5RA > I for some A G a{JF)- 
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The finite difference metfiod is calied elementary stable if, for any value of the time step At, the linear 
stability of each ^„ G is the same as the stability of each ^ S F^. 

From a theoretical point of view this criterion is good in every case to control the stability of the 
steady-states. Nevertheless, in the case of multi-dimensional systems, expecially in two and three dimensions, 
it is more convenient to use an alternative method, the so-called Jury criterion Let Jp •^F(Cn)i 
two dimensional systems the characteristic polynomial of the Jacobian can be written as: 

- XtrJp + detJp (5) 

in order to have iRA < 1 for all A G a{JF), the Jury condition states that: 

\trJF\<l + detJF<2. (6) 

Therefore, this criterion establishes that exists a necessary and sufficient condition to guarantee the asymp- 
totical stability of the steady-state solutions. 

The forward-Euler method. This is one of the oldest way to derive a finite difference equation from a 
differential equation. This numerical procedure requires that the system (1) is transformed by introducing 
these substitutions: 

m (7) 



d(,{t) — ^„ 

^ ^ At ' 

where At is the step size and £,n+i ~ £,{t + nAt). As showed by Mickens (see i.e. the reference [ ' in many 
cases the transformation of eq. (8) leads to numerical instabilities which occur as solutions of the finite 
difference equations, steady-states, bifurcations etc. that do not appear in the originary differential equation. 

3.1 Non Standard Finite Difference method 

To avoid these potential problems, Mickens [ t_] suggests what is known as the Non Standard Finite Difference 
(NSFD) method. An important concept on which the NSFD schemes is based is the Dynamic Consistency 
[17]. This criterion states that, let the first-order autonomous ODE 

f^/(e,fc) (9) 

and given U its set of properties, the difference equation 

e„+i =F(C„,fc,At) (10) 

is dynamically consistent with (9) if it respects the same set of properties. As indicated in [I I], the dynamic 
consistency is verified when a difference equation possess the same stability, bifurcations and chaos of the 
original differential equation. 

A finite difference method is called a NSFD scheme is almost one of the following conditions is satisfied: 

L The order of the discrete derivative should be equal to the order of the corresponding derivatives of the 
differential equation. 

IL Denominator functions for the discrete rapresentation must be non-trivial. The following replacement 
is then required: 

At — > (j){At) + 0{At^) 
where 0(Ai) is such that < 0(Ai) < 1 MAt > 0. 
III. Nonlinear terms must be replaced by non-local discrete representations, i.e. 
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IV. Special conditons that hold for the solutions of the differential equations should also hold for the 
solutions of the finite difference scheme. 

V. The scheme should not introduces spurious solutions. 

The NSFD scheme incorporates the principle of Dynamical Consistency. 

A very important characteristic of dynamical systems, especially in those of biological interest, is that 
all solutions must remain non-negative in order to maintain the problem well-posed, from biological and 
mathematical points of view. A method that respects the Mickens rules (I to V) and preserves the solution 
positivity is called Positive and Elementary Stable Non standard (PESN) method. The appropriate non-local 
approximation that provides a such scheme is constructed by using the indications suggested by Patankar 
[20]. Applying these measures, the discrete scheme will comply with the physical properties of the differential 
equations, without any restriction on the step size At. 



4 An elementary application: the logistic population model 

The logistic equation, originally due to Pierre Francois Verhulst [ ], is the law that regulates with good 
approximation the growth rate of a certain population number as function of time. The model is based on 
the assumption that the population evolves in an environment with limited resources with no immigration or 
emigration phenomena. Let x{t) the population at instant t, the law that regulates it can be expressed from 
the following first-order autonomous ODE: 

= rx{t) - ^ (11) 
dt ^ ' E ^ ^ 

where E is the carrying capacity of the system and r is the intrinsic growth rate {r = h — d, where h is the 
instantaneous birth rate and d the instantaneous death rate); also we assume r, fc > 0. The eq. (11) is the 
same as (1) assuming ^ — x{t) and / coinciding with the right hand term. In this equation the term rx{t) is 
the intrinsic population growth rate, while —rx'^{t)/E is the self- limit at ing term. These two factors of / have 
opposite signs then / can be greater, lesser or equal to zero: this determines respectively the increase, the 
decrease or the constancy of the population. 

Applying the condition dieted by eq. (2), is simple to calculate the two steady-states: a trivial solution 
ii = and a non-trivial X2 = E. Sobstituting these solutions in the Jacobian J{x) = \ — 2x/E the following 
inequalities J(ii) > and J{x2) < hold. We can deduce that only the solution X2 — E is asimptotically 
stable. 

This can be deduced also by an alternative way. The differential equation (11), solved by variables separation 
technique, admits the following solution: 

xit) = , s -■ (12) 

Calculating the limit of eq. (12) for t — > +oo the population has just E as its asymptotic value. 
4.1 Discrete model of logistic equation. 

Following the NSFD rules, the logistic model of eq. (11) can be transformed in the discrete scheme: 

^-rx(t)-'-^ (13) 

where Ax(t) — — Xn and '/'(At) is a function of the time step in respect ofthe rule II; in accord with 
rule III, — > XnXn+i- Therefore Xn+i can be so explicited: 

l + r4>{At) 
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where xq = x{t = 0) is the initial population. It is simple to show that the condition of positiveness Xn+i > 
is automatically satisfied if cj) satisfy the rule II of the NSFD scheme. The steady-states of eq. (14) are Xi — 
and xi = E: for the trivial steady-state, J(ii) = 1 -I- r0 > 1 for r > and < (/) < 1, then is unstable; for the 
non trivial one, J{x2) — 1/(1 -I- re/)) < 1, then is asimptotically stable. We conclude that the NSFD scheme 
well reproduces the behavior of the differential equation (13). 

In opposition, the Euler and Runge-Kutta discretization methods are dynamically inconsistent with the 
original differential equation [] >!], showing numerical instabilities. Also in [18] is remarked that the application 
of these two standard schemes produces a not correct monotonicity and oscillations for particular values of 
the mesh At. 



5 Mutualistic interaction between two species 

It is now considered a more complex case respect to the logistic equation: this is the formulation introduced 
by Holland and DeAngelis 2009 [:!] for the mutualism description, based on Rosenzweig-MacArthur model 
of predator-prey interactions. After a brief review of the problem in continuous-time domain, two discrete 
schemes will be proposed. 

5.1 Continuous-time domain 

Subsequently it will be denoted with x{t) and y{t) the numbers (or densities) of the mutualistic species. In a 
mutualistic model, the growth rate of each species must involve the same terms of the logistic equation (11) 
that represent effect separate from mutualism: a linear term that represent the intrinsic population growth 
rate and a quadratic one to modify the growth rate with density dependent self-limitation. Also, for each 
species, in case of bi-directional mutualism, other two terms are need: 

• a positive term that quantify the advantages that the presence of a species induces on the growth of 
the other (aij in eq. (15)), 

• a negative term that describes the disadvantages that the presence of a species induces on the growth 
of the other {f3i in eq. (15)). 

With these two conditions, a community of two mutualistic species can be described by the following system 
[3]: 



dy 
It 



r2y + a2i^~P2^ 



2+y 



(15) 



I x{0) > 0, y(0) > 



where all parameters are constants, reals and non-negative. To create correspondence with the notation of eq. 
(1), we assume ki — (ri, ai2, Pi, ei, di) and ^2 = (r2, a2i, hi, /32, 62, ^2). 



5.2 Discrete-time domain 

Forward Euler method. It is required to apply the rules (7) and (8) to the differential equation system 
(15). The result is: 



At ~ 



riXn + ai2i 



h2+yn 

-'^22/„ + a2l7f^ 
, xo > 0, > 0. 



^2 e2+y„ 



- dixl 
d2yl 



(16) 
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that can be so explicited: 



(l + At (n - + 2/„ - ^))) 

Vn+l ^yn{l + A< (r2 - d^Vn + X„ (^f^ - 



(17) 



xo > 0,?;o > 0. 



It is worth noting that the above system is not unconditionaUy positive due to the presence of negative terms 
that can allow the generation of negative solutions. This feature is considered to be precursor for numerical 
instabilities. 

For the calculation of the fixed points, the following system must have resolved: 



Vn = Vn 



_§2 



(18) 



, 2^0 > 0, 2/0 > 0. 



The scheme (17) and the solutions coming from system (18) will be compared in the Section C with results 
obtained by using NSFD scheme that will be shown in the following paragraph. 

NSFD scheme. Now, we transform the differential equations in system (15) respecting the following NSFD 
rules: 

• II rule: denominators of left-hand side are replaced with (j) = (j){At); 

• III rule and PESN criterion: the terms on the right of the equations are replaced with a non-local 
approximation. 

Then, the resulting finite difference system is: 

riXn + Oli2-^ 



d2ynV 



n+1 



(19) 



. a;o > 0, yo > 0. 
Then they are calculated Xn+i and ijn+i- 

Xn+l = -Fl(x„,?/„,0, fci) 

, xo>o,yo> 0. 



(20) 



where: 



F2{Xn,yn,(l),k2) 



(ei + Xn){h2 + Vn) + '/'(ei + Xn){ri{h2 + Vn) + ;/««i2) 
(ei + a;„)(/i2 + Vn) + (t>{h2 + Vn){dlXn{ei + Xn) + VnPl) 

(62 + yn){hi + Xn) + (t>{e2 + yn){r2{hl + X„) + X„0;2l) 
(62 + yn){hl + Xn) + 4>[hi + Xn){d2Vn{e2 + Vn) + Xnl32) 



(21) 
(22) 
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in which fci and k2 are the parameters already defined in Subsection 5.1. Since all components of ki and k2 
are non-negative we have Xn+i,yn+i ^ for all n > 0, guaranteeing the respect of the positiveness condition. 
Introducing f2 ~ lj|gp^ (t( J), (j) is given by the following expression [1-1]: 

0(At) = '- (23) 

q 

in which: 

Following the definition given in (4), the steady-states are calculated by solving the below system: 

Fl{Xn,yn,(l),kl) = Xn 

(25) 

F2{Xn,yn,(t>,k2) = Vn 

To discuss the stability of all ^„ G Td^ it is convenient to check if the Jury condition mentioned in eq. (6) 
it is respected for each of these. To do this, it is need to calculate the Jacobian matrix Jpixmyn) — (jjj)2x2- 
The elements of matrix are the following: 

, y ) - ^^^ (f'ri)ih2 + y„) + (?!)2/„ai2((ei -I- x„)'^ + (j)ynPi{ei + 2x„)) 

. , _ 0a:„(ei -I- a;„)(/i2(ei -I- a:„)(l-|- 0dix„)ai2 - ((1 + 0ri)(/i2 -I- j/„)^ -I- (^y^ai2)/3i) 

^"^ " {h2 + 2/„)2((ei + x^){l + M:e„) + 0y„/3i)2 ^^^^ 

^ ^ _ '/'yn(e2 + ^,0(^*1(62 + yn)(l + </'fi2yn)Q;21 - ((1 + </'»^2)(/ll + a;„)^ -I- 0X^Q;2i)/32) 

(/ii + a;„)2((e2 -I- 2/„)(l-l- (/)(i2yn) + 02:„/32)2 
, y )-''"'" ^ (t'r2){hi + Xn) + 0a;„Q;2i((e2 + ?/„)^ + 0.T„/32(e2 + 2y„)) 

6 Numerical simulations 

In this Section are presented some numerical simulations to test the dynamic consistence of the NSFD and 
the forward Euler schemes found in Subsection 5.2 with the continuos-time one described by the system (15). 

6.1 Steady-states analysis 

In continuous-time domain the problem is largely discussed in [1] and [ ]: these authors have performed 
many simulations assuming different values for the model's parameters, i.e. assuming different types of 
mutualism. In this paper a unique simulation for steady-state analysis is reported, hypothesizing for the 
involved parameters ri, r2, ai2, Q!2i, 132, ei, 62, di and c?2 the same values used by Holland and DeAngelis 
[.-',] for bi-directional mutualism. The steady-states have been calculated applying the condition (2) to the 
system of differential equation (15). From graphical point of view, the steady-states are found plotting the 
zero-growth isoclines of each system's equation and are represented by the points where the isoclines intersect 
between them. The phase- plane diagrams for the x-y dynamics is showed in Figure la where the zero-growth 
isoclines are depicted. 

In the discrete case, NSFD derived scheme has been compared with the Euler one calculating the 
solutions of the systems (25) and (18) respectively. The Figure lb represents the isoclines described by the 
Euler and NSFD schemes for At = 10~^. From those graphs can be seen that, for the assumed parameters, 
the dynamical system has 6 steady-states, indicated with F", F^,.., F^ in counterclockwise sense. Also, none 
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Figure 1: Comparison between phase-plane diagrams for continuos-time and discrete-time cases, (a) Zero-growth isoclines for 
system of eq. (15) are represented as black {dx/dt = 0) and red {dy/dt = 0) lines, (b) Zero growth isoclines come from the 
system (25): black dotted line for first equation and red dotted line for the second equation; zero-growth isoclines come from 
system (18): in blue the first equation and in green the second one (completely masked by NSFD lines). For the discretization 
parameters, the following values have been used: At = 10~^, q = 10. The steady-states are marked with letters T'^ to F^: the 
square (F", F^, F'*) indicates the repellers of the system as well as the circle (F^, F'', F^) indicates the attractors. The value 
assigned at each parameter is [ >]: ri = r2 = 0.3, «i2 = «2l = 0.6, qi = q2 = 1, P\ = h = 0.2, di = d2 = 0.01, /ii = /i2 = 0.3 
and ei = 62 = 0.3. 



Steady-state 


x-coordinate 


y-coordinate 


Ai 


A2 


Stability 


pO 








0.30 


0.30 


R 


F^ 


30 





-19.11 


-0.30 


A 


F^ 


25.62 


83.17 


-0.80 


0.41 


R 


r3 


69.83 


69.83 


-0.70 


-0.30 


A 


r4 


83.17 


25.62 


-0.80 


0.41 


R 


rb 





30 


-19.11 


-0.30 


A 



Table 1: Steady-states of Figure 1. A=attractor, R=repeller. 



detJ + 1-\trJ\ detJ-1 




Figure 2: (a) Plot of eq. (27) versus the time step At. Each curve is related to at least one of the steady-states F", F^, F^. 
(b): plot of eq. (28) versus At. In both figures, the solid curves are related to the NSFD and the dashed curves to the Euler 
method. 
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of the two discretization methods introduces spurious solutions: the coordinates of the steady-states are the 
same and the isochned are perfectly superimposed. The Table 1 summarizes the results of these simulations 
and the conclusions on steady-states stability for continuous and discrete-time schemes. From the examination 
of the signs of the eigenvalues, it emerges that r°, and F^ are repellers (indicated in figure with a black 
square) and F^, F^ and F^ are attractors (black circle). 

However, when At is increased up to 8.445 • 10~^ we observe that the Euler scheme converges toward 
wrong steady-states or the two populations have a boundless growth. On the contrary the phase-plan of the 
NSFD remains unchanged. 

Since one of the requisite for the dynamic consistence is the invariance of the stability conditions regardless 
to the discretization step-size, has been checked for NSFD scheme if the Jury condition is guaranteed or not 
independently from the value of At. This check is not performed for Euler scheme because, as explained 
above, it fails the convergence to the correct fixed points for certain values of the step-size. 

In order to do this verify, the condition (6) is splitted in two inequalities what will be proved for all 
5 S Fd, rendering explicit the dependence from At of the Jacobian elements: 



The Figures 2a and 2b represent respectively the left term of eq. (27) and (28), both evaluated as function of 
At: in both graphics every curve is related to at least one steady-state F°, F^, ...,F^. Only for reasons of 
correct graphic visualization the two figures have different scales on abscissa. It is easy to see that both the 
inequalities (27) and (28) are satisfied for NSFD schemes only from F^, F^ and F^ steady-states for each 
value of At represented. Nevertheless, has been checked that all the curves are monotone for < Ai < 10, 
then the conclusions about the steady-states stability are unchanged up to At = 10. 

6.2 Transient dynamics 

In this Section a transient analysis has been performed with the purpose to emphasize the possible criticalities 
of the two considered discretization schemes. 

The phase-plan portraits of NSFD and forward Euler difference equations has been showed in Figure 3 a 
and b respectively: all the used parameters are the same used for computation of Figure 1, exception for the 
time step that now is At = 6 • 10^^. The red points are the values assumed for the initial populations and 
blue points represent the evolution of the system with the increase of iterations. From the comparison it is 
possible to observe some differences highlighted by red rectangles in the Figure 36 named with Di and D2. 
The Figures 3c and d represent the magnification of what happens in Di and D2 respectively: it is possible 
to note that in subdomain Di the population x assumes negative values for about y > 45 and the population 
y assumes negative values for about x > 45. All these negative values appear in the transient. 

The Figures 4a and c show a plot of the population vs the iteration number using respectively (17) and 
(20) schemes. In order to highlight the numerical instabilities, larger values of At are used in this computation: 
respectively 8 • 10~^ for 4a and 8.44 • 10^^ for 4c; the initial population size is (xo,yo) = (30, 10). It can be 
seen from the Figures 4a that for NSFD equations, after a transient in which the populations x and y reach 
their maximum value about 72 and 11 individuals, they converge smoothly toward 30 and individuals. For 
Euler equations there are at least two undesiderable features: the jump after 40 iterations about and some 
negative values after about 60 iterations, as showed in zoomed window (green line). Nevertheless, with both 
methods, all the populations converge to the correct steady-state F^. The relative phase-plan is depicted in 
the Figure 4&, where the different trajectories are evident, especially in the first iterations. 

For computation of Figure 4c the conditions are the same, exception for the time step that now is 
At = 8.44 • 10~^. The remarkable fact is that, while for NSFD the transient preserves the same aspect of 
Figure 4a, the Euler scheme now converges toward F^ instead of F^. In Figure id is plotted the relative 
phase-plan where the anomaly is clear. 

Is important to note that none of these anomalies occur with At < 10~^: in particular can be easily 
checked that with At = 10"'^ the curves of evolution of the populations coming from NSFD and Euler schemes 
are perfectly superimposable. The first numerical instabilities in forward Euler scheme begin to occur when 
At overcomes the critical value of Ate ~ 5 • 10~^. The reason is that in forward Euler, the step sizes larger 



detJ{At) - 1< 
detJ{At) + 1 - |trJ(At)| > 



(27) 
(28) 
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Fi gure 3: (a) Phase-plane portrait of system (20). (b) Phase-plane portrait of system (17). (c) Particular of subdomain 
Di belongs to (b). (d) Particular of subdomain D2 belongs to (b). Has been used the same parameters used for Figure 1 
computation, At = 6 ■ 10~^. 




100 150 200 

iterations 



30 40 
X 



(d) 



Figtire 4: Comparison between transients of Euler, eqs. (17), and NSFD difference equations, eqs. (20). (a) x populations 
vs iteration number: red points for NSFD and blue points for Euler scheme; y populations vs iteration number: violet points 
for NSFD and green points for Euler scheme. At = 6 ■ 10~^. (b) Phase-plan of (a), (c) The same of plot (a), but with 
At = 8.44 ■ 10^^. (d) Phase-plan for (c) simulation. 
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than Afc are greater than some characteristic time relevant of this physical problem. The same is not true 
employing the NSFD technique, in which the numerical instabilities are eliminated since the selected step 
size of eq. (23) is never larger than the smallest time scale of the system. 

7 Conclusions 

On the base of the work published by Holland and DeAngelis in 2009, in which is constructed and analyzed a 
population dynamic continuos-time model that link consumer functional responses of one mutualistic species 
with resources supplied by another, in this paper a strong non-standard finite difference scheme to solve the 
continuous-time model is constructed. As showed by numerical simulations, the NSFD framework reproduces 
the dynamical features of the continuos-time analogous. In particular, it is proved that the method guarantees 
the correct asymptotic behavior regardless from the size of the time step. It was shown that using the same 
initial conditions the forward Euler induces scheme-dependent numerical instabilities for At greater than a 
critical values, that manifest as jumps in the population evolution in time domain and convergence toward 
wrong steady-state solutions. On the other hand, in tests that were performed, the NSFD scheme shows none 
of these problems confirming the versatility of the method. 
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